library(plyr)
library(estimatr)
library(emmeans)
library(tidyverse)
library(magrittr)
library(broom)

d1 <- "https://github.com/thomasjwood/full_fact/raw/main/ff_mt.rds" %>% 
  url %>% 
  gzcon %>% 
  readRDS %>% 
  mutate(
    out_agree = out_agree %>% 
      mapvalues(
        c("Tend to disagree",
          "Tend to agree"),
        c("Disagree",
          "Agree")
      ) %>% 
      factor(
        c("Strongly disagree",
          "Disagree",
          "Neither agree nor disagree",
          "Agree",
          "Strongly agree")
      ),
    out_true = out_true %>% 
      factor(
        c("True",
          "Probably true",
          "Not sure",
          "Probably false",
          "False") %>% 
          rev
      ),
    agree_num = out_agree %>% 
      as.numeric,
    true_num = out_true %>% 
      as.numeric,
    cond = cond %>% 
      factor(
        c("cond_items",
          "cond_misinfo",
          "cond_corr")
        )
    )

d1$com_num <- d1$agree_num %>% 
  add(
    d1$true_num
  ) %>% 
  divide_by(2)

d1$ideo_num <- d1$ideo %>% 
  mapvalues(
    c("Left1", 
      "1",
      "2",
      "3",
      "4",
      "5",
      "6",
      "8",
      "7",
      "9",
      "10",
      "Right10",
      "Prefer not to say", 
      "Don't know",
      "DK"),
    c(1, 1:10, 10, NA, NA, NA)
  ) %>% 
  as.character %>% 
  as.numeric


# setting up modelling table ----------------------------------------------

t1 <- d1 %>% 
  mutate(
    country = issue %>%
      str_detect(
        "global|saltwater"
      ) %>%
      ifelse(
        "Multi-country",
        country
      )
    ) %>% 
  filter(
    wave == "wave 1"
  ) %>% 
  group_by(
    country, issue
  ) %>% 
  nest

t1$mods <- t1$data %>% 
  map(
    \(i)
      
      lm(
        com_num ~ cond * ideo_num,
        data = i
      )
    )

t1$emm <- t1$mods %>% 
  map(
    \(i)
      
      i %>% 
      emmeans(
        consec ~ cond | ideo_num, 
        at = list(
          ideo_num = seq(1, 10, length.out = 500)
        ),
        level = .95
      )
  )


t2 <- t1 %>% 
  ungroup %>% 
  select(
    country, issue, emm,
  ) %>% 
  pmap_dfr(
    function(
      country, issue, emm
    ){
      
      j <- emm %>% 
        extract2("contrasts") %>% 
        tidy %>%
        as_tibble %>% 
        filter(
          contrast == "cond_corr - cond_misinfo"
        ) %>% 
        mutate(
          sig = adj.p.value %>%
            gtools::stars.pval() %>% 
            as.character
        )
      
      rlj <- j$sig %>% 
        rle
      
      j$sg <- "g_" %>% 
        str_c(
          rlj$values %>% 
            seq
        ) %>% 
        rep(
          rlj$lengths
        )
      
      j %>% 
        mutate(
          issue = issue,
          country = country
        )
    }
  )

t2 %<>% 
  mutate(
  ymin = estimate %>% 
    subtract(
      std.error %>% 
        multiply_by(1.96)
    ),
  ymax = estimate %>% 
    add(
      std.error %>% 
        multiply_by(1.96)
    )
  )

t2$issue %<>%
  factor(
    t2 %>% 
      group_by(
        issue
      ) %>% 
      nest %>% 
      mutate(
        estimate = data %>% 
          map_dbl(
            \(i)
              lm(
                estimate ~ ideo_num, data = i
                ) %>%
              tidy %>% 
              filter(
                term == "ideo_num"
              ) %>% 
              use_series(estimate)
            )
        ) %>% 
      arrange(desc(estimate)) %>% 
      use_series(issue) %>% 
      as.character
    )


t2$lab <- t2$estimate  %>%  
  round(1) 

t2$lab <- t2$lab %>% 
  str_detect("\\.") %>% 
  ifelse(
    t2$lab,
    t2$lab %>% 
      str_c(
        ".0"
      )
  )

t2$lab %<>% 
  str_c(
    gtools::stars.pval(t2$adj.p.value)
  )

t2$sig %<>% 
  mapvalues(
    c("*", "**", "***"),
    c("* p < .05",
      "** p < .01",
      "*** p <.001")
    ) %>% 
  factor(
    c("* p < .05",
      "** p < .01",
      "*** p <.001")
    )

t2$issue %<>% 
  mapvalues(
  c("nigeria_youthue", "nigeria_birthreg", "safrica_tobacco", "nigeria_chancellor", "nigeria_malaria",
    
    "haigh_knife", "pichetto_shanty", "safrica_latrines", "ashworth_nurses", "fernandez_debt",
    
    "fernandez_killings", "safrica_civserv", "nigeria_crossriver", "macri_debt", "macri_arrests",
    
    "safrica_resbank", "cele_crimestats", "johnson_educ", "global_cooling", "vaccine_cancer",
    
    "javid_homeless", "saltwater_covid"),
  
  c("Nigeria-Youth U/E 70%", 
    "Nigeria-70% births not registered", 
    "S.Africa-49% men, 34% women, smoke tobacco",
    "Nigeria-Nigerian unis 15 female chancellors",
    "Niveria-Malaria kills 300k annually",
    
    "UK-Knife crime all-time high",
    "Argentina-4k shanty towns in Buenos Aires",
    "S.Africa-100s children drown in pit latrines",
    "UK-200k nurses quit NHS since 2010",
    "Argentina-Debt increased 38% to 100% GDP",
    
    "Argentina-3,262 murders last year",
    "S.Africa-80% of govt spending is salaries",
    "Nigeria-Cross River lowest national crime",
    "Argentina-Debt is today 100% GDP",
    "Argentina-85k arrested for drug trafficking",
    
    "S.Africa-German owns 57% SA Reserve Bank",
    "S.Africa-SA only country releases crime stats",
    "UK-Johnson invests £14bn in education",
    "Multi Country-Global cooling '16-18",
    "UK-Vaccines might cause cancer",
    
    "UK-Homeless halved since 2008",
    "Multi Country-Saltwater kills COVID"
  )
)

t2 %>% 
  ggplot() +
  geom_ribbon(
    aes(
      ideo_num, ymin = ymin, ymax = ymax
    ),
    size = .25,
    color = "grey1",
    fill = "grey95"
  ) +
  geom_ribbon(
    aes(
      ideo_num,
      ymin = ymin,
      ymax = ymax,
      fill = sig,
      group = sg
    ),
    size = .25,
    data = t2 %>%
      filter(
        adj.p.value <= .05
      ),
    color = "black"
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dotted"
  ) +
  geom_point(
    aes(
      ideo_num, estimate
    ),
    data = t2 %>% 
      group_by(
        issue
      ) %>% 
      slice(
        c(1, 250, 500)
      ),
    shape = 21,
    size = 8,
    color = "grey5",
    fill = "grey95"
  ) +
  geom_text(
    aes(
      ideo_num, estimate, label = lab
    ),
    data = t2 %>% 
      group_by(
        issue
      ) %>% 
      slice(
        c(1, 250, 500)
      ),
    size = 1.8, 
    family = "Roboto"
  ) +
  facet_wrap(
    ~issue, 
    labeller = label_wrap_gen(width = 20),
    ncol = 6
    ) +
  scale_fill_grey(
    start = .8, end = .4, 
    guide = guide_legend(
      direction = "vertical",
      title.position = "top"
    )
  ) + 
  scale_x_continuous(
    breaks = c(
      1.25, 5.5, 9.75
    ),
    labels = c("Left",
               "Mod.",
               "Right"),
    expand = c(.1, .1)
  ) +
  scale_y_continuous(
    expand = c(.05, .05),
    breaks = seq(
      -1.5, .5, .5
    ),
    labels = seq(
      -1.5, .5, .5
      ) %>% 
      str_replace_all("0\\.", "\\.")
    ) +
  labs(
    x = "Ideology",
    y = "Correction effect",
    fill = "Effect significance"
  ) +
  theme(
    strip.text.y = element_text(angle = 0),
    legend.margin = margin(-.1, 0, 0, 0, "cm"),
    panel.grid.minor = element_blank(),
    plot.caption = element_text(
      face = "italic", size = 7,
      ),
    legend.position = c(.875, .125)
  )
